An analytic model for the bispectrum of galaxies in redshift space 
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We develop an analytic theory for the redshift space bispectrum of dark matter, haloes and 
galaxies. This is done within the context of the halo model of structure formation, as this allows for 
the self-consistent inclusion of linear and non-linear redshift space distortions and also for the non- 
linearity of the halo bias. The model is applicable over a wide range of scales: on the largest scales 
the predictions reduce to those of the standard perturbation theory (PT); on smaller scales they 
are determined primarily by the nonlinear virial velocities of galaxies within haloes, and this gives 
rise to the [/-shaped anisotropy in the reduced bispectrum - a finger print of the Finger-Of-God 
distortions. We then confront the predictions with measurements of the redshift space bispectrum 
of dark matter from an ensemble of numerical simulations. On very large scales, k = 0.05 frMpc -1 , 
we find reasonably good agreement between our Halo Model, PT and the data, to within the errors. 
On smaller scales, k — 0.1/iMpc _1 , the measured bispectra differ from the PT at the level of 
~ 10 — 20%, especially for colinear triangle configurations. The Halo Model predictions improve 
over PT, but are accurate to no better than 10%. On smaller scales k — 0.5 — 1.0 /tMpc -1 , our model 
provides a significant improvement over PT, which breaks down. This implies that studies which 
use the lowest order PT to extract galaxy bias information are not robust on scales k > 0.1 feMpc -1 . 
The analytic and simulation results also indicate that there is no observable scale for which the 
configuration dependence of the reduced bispectrum is constant — hierarchical models for the higher 
order correlation functions in redshift space are unlikely to be useful. It is hoped that our model 
will facilitate extraction of information from large-scale structure surveys of the Universe, because 
different galaxy populations are naturally included into our description. 

PACS numbers: 98.80.-k 



I. INTRODUCTION 



Statistical analyses of the large scale structures ob- 
served in galaxy surveys can provide a wealth of in- 
formation about the cosmological parameters, the un- 
derlying mass distribution and the initial conditions of 
the Universefl 1 1 1 1- Some complex combina- 
tion of the information is commonly extracted through 
measurement of the 2-point correlation function, or it's 
Fourier space analogue the power spectrum. Since the 
evolved density field of galaxies is highly non-Gaussian, 
further complementary information is contained within 
the higher order clustering statistics For example, 

analysis of the large-scale 3-point correlation function, 
or its Fourier space dual the bispectrum, on large scales, 
can: test the non-linearity of bias - the way in which an 
observable tracer distribution samples the unobservable 
distribution of physical interest [1 H El EH ; constrain 
our hypothesis of Gaussianity in the initial conditions 

EE EI El; break de 

generacies between parameters, 
hence allowing improved constraints on the amplitude of 
the matter power spectrum. In addition, higher order 
statistics have been highlighted as an important piece of 
solving the puzzle as to whether the observed accelerated 
expansion of the Universe is due to Dark Energy physics 
or a modification to gravity [H E3j EE EH ■ On smaller 
scales these statistics can be most usefully used as a dis- 



criminator for the shapes of haloes [20( - and thus have 
the potential to constrain the small-scale dark matter 
physics. 

The current state of the art galaxy redshift surveys 
[2ll I22I I23I |24| have provided large samples of the Uni- 
verse, and investigators have already carried out some of 
the tests noted above: 3-point correlation functions have 
been estimateclby [11, EE E3, El El EH and bis P ectra 

by m muni- 

The recovery of precise cosmological information from 
these measurements is not straightforward, owing to the 
influence of non-linear mass evolution, biasing and red- 
shift space distortions. Current theoretical modeling of 
the large scale bispectrum rests on results from pertur- 
bation theory and non-linear biasing in real (as opposed 
to redshift) space. In the large scale limit, this gives rise 
to a simple model d [1 EH : 

Q s (ki,k 2 ,9 12 ) = l[Q m (fc 1 ,A: 2 ,0 12 ) + C f] ; (1) 

where the functions Q l are the reduced bispectra of galax- 
ies and matter, Qj 23 = B l 123 / [PfP^ + P'P* + PiP[} . In 
the above -B123 = B{k%, k 2 , #12) is the matter bispectrum 
and P a = P{k a ) the matter power spectrum. The co- 
efficient 6f is the large-scale linear bias parameter and 
c| = h\jb\ is the first non-linear bias parameter. It is 
usually assumed that this relation holds in redshift space 
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as well, but that it does not in detail can be seen from the 
work of [3f| . That it nevertheless appears to be a reason- 
able working hypothesis was demonstrated by [36| , We 
shall, hereafter, refer to this model in real and redshift 
space as PT and PTs, respectively. 

There are several reasons why we wish to improve 
upon PTs. Firstly, it is well known that in redshift 
space the distortion effects from non-linear structures, 
such as Finger-of-God (hereafter FOG) distortions, pol- 
lute the scales that are usually identified for linear treat- 
ment. This can not be accounted for in the perturbation 
theory in any other way than supposing 'ad-hoc' fixes 
to the model [H, [37|]. A pragmatist might argue that 
one may take scales that are sufficiently large that these 
corrections can be neglected. However, even if we are 
proficient enough to accurately separate linear from non- 
linear scales, then we are still faced with loosing a signif- 
icant amount of information from our data through the 
restrictions to very large scales. Therefore some means 
for robustly modeling the FOG effects is clearly of great 
value as this may allow us to expand the utility of our 
data set and improve precision. 

Secondl y, in our study of the large-scale galaxy power 
spectrum [381 ] we found that there was non-trivial scale 
dependence arising from non-linear bias and gravita- 
tional mode coupling, even on the largest scales currently 
probed. One may then ask how these properties affect 
the predicted bispectra. 

Thirdly, if we assume that that galaxy velocity field is 
an unbiased tracer for the velocity field of dark matter, 
then through studying the higher order clustering statis- 
tics in redshift space we have a direct probe of the sta- 
tistical information of the dynamics of the CDM density 

field itself [Si, SO, El,!!!. 

In this paper we build a new analytic model for the 
fully non-linear redshift space bispectrum. We will con- 
centrate on the isotropically averaged (i.e. the monopole) 
bispectrum, for reasons of simplicity. We work in the con- 
text of the halo model [43j, since it naturally affords a 
means for includin g lin ear and non-linear density and ve- 
locity information [3, SI, liE lijl. liil and neatly allows 
for the inclusion of galaxies "EiniallO, HH, [H, HI). Fur- 
thermore, as was shown in [38| the halo model presents 
a natural framework for understanding the origins of the 
non-linear scale dependence of bias. However, the limita- 
tions of the halo model predictions for precision measure- 
ments of the matter power spectrum on lar ge sc ales have 
been known for some time now [HI, H3, 155ll5t| . We will 
therefore use measurements from numerical simulations, 
to confirm the validity of our predictions. 

The paper breaks up as follows: In Section [H] we for- 
malize the halo model in redshift space, providing general 
expressions for the 3-point function and bispectrum. Sec- 
tion IIIII details the necessary components of the model; 
we pay special attention to the redshift-space clustering 
of halo centers. Section IIVI presents the central ana- 
lytic result of the paper — a calculation of the bispectrum 
monopole. Some results of evaluating our expressions are 



presented in Section |V1 In Section I VII we confront our 
model with measurements from iV-body simulations. In 
Section [VIII we summarize our conclusions. 

Although our analysis is general, we shall illustrate 
our results with specific examples. When necessary, 
we assume a flat Friedmann-Lemaitre- Robertson- Walker 
(FLRW) cosmological model with energy density at late 
times dominated by a cosmological constant (A) and a sea 
of collisionless cold dark matter particles as the dominant 
mass density. We set fi m = 0.27 and Qa = 0.73, where 
these are the ratios of the energy density in matter and 
a cosmological constant to the critical density, respec- 
tively. We use a linear theory power spectrum generated 
from cmbfastJH^, with baryon content of = 0.046 
and h = 0.72. The normalization of fluctuations is set 
through (78 = 0.9, which is the r.m.s. variance of fluctu- 
ations in spheres of radius 8/i -1 Mpc. 



II. HALO MODEL IN REDSHIFT SPACE 

A. Formalism 

In the halo model (see [43[ for a review) the density 
field is decomposed into a set of dark matter haloes, 
where a halo is defined to be a region that has under- 
gone gravitational collapse forming a dense virialized ball 
of cold dark matter (CDM). All statistical quantities of 
interest are then considered as sums over the halo dis- 
tribution. Thus to understand the large scale clustering 
of a distribution of objects, haloes, galaxies or dark mat- 
ter, we simply require understanding of how the haloes 
themselves cluster; the different tracer types simply act 
as weights. In particular, different galaxy populations 
'weight' haloes differently: the Halo Occupation Distri- 
bution (HOD) [US Hi HHH, specifies how the prob- 
ability for obtaining N galaxies depends on halo mass M. 
To model redshift space statistics, we require additional 
information about how the large scale velocity field mod- 
ifies the halo clustering, as well as a model for the dis- 
tribution function of galaxy velocities within each halo. 
The halo model in redshift space, at the 2-point level, was 
developed by [IE SE S3 (hereafter we shall refer to the 
Halo model in real and redshift space as HM and HMs, 
respectively). However, some unresolved issues remained 
with regard to the base formalism. These were resolved 
by [48l ] and our description of 3-point statistics presented 
here extends these analyses. For completeness some of 
these details are repeated below. Before continuing, we 
note that the problem of redshift space distortions in the 
Halo Model was also recently addressed by |58j , who used 
numerical simulations to construct an empirical model 
for the distribution function of halo pair- wise velocities. 
Our approach is complimentary to that, since the results 
for the large-scale halo clustering are derived within the 
context of the analytic perturbation theory rather than 
being fit for. 

The density field of dark matter, haloes or galaxies 
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may be written as 

p s a (s)=J2[Wa] i K ji ( S - Si \M i ) , (2) 

i 

where a = {1,2,3} refers to the particular choice of 
weight for the ith halo in the sum, i.e. [W a ]i — 
{1, iV g (Mj), . . . } depending on the spectra one 
wishes to model, and where N 6 (Mi) is the number of 
galaxies in halo i. i is the normalized density distri- 
bution of objects in redshift space within the ith halo. In 
this paper we shall always assume that t — p s (s)/M, 
is the mass-normalized density profile of dark matter in 
redshift space, although our formalism does not rely upon 
this assumption and may readily be generalized for more 
complicated mass distributions. At this point the only 
difference between Eq. ^ and the real space density 
field, is that we have used s to denote comoving spa- 
tial positions. However this has the special meaning that 
Hubble's law, v = H(a)r, is used to infer proper ra- 
dial positions from recession velocities, where v is the 
proper velocity, H(a) = a/ a is the Hubble parameter 
and r is the proper separation (related to comoving co- 
ordinate through r = ax). The notion of redshift space 
distortions then follow from the fact that objects which 
form through gravitational instability acquire a local pe- 
culiar velocity of their own, and hence the velocity-space 
mapping in general is non-linear. In this paper we shall 
work in the plane parallel approximation, where observed 
structures are located at infinity. Then the mapping is 

s z = z - u z (x) ; s_l = x_l , (3) 

I 



Cq(si,s 2 ,s 3 ) 

Cih( S 1, S 2,S 3 ) 
Ca,2H( S l> S 2,S 3 ) 



Ca,3H( s l> s 2,S3) 



where ££ c and (f[ c are the 2- and 3-point correlation func- 
tions of halo centers, conditioned on halo masses and 
where n(M)dM is the halo mass function, which gives 
the number density of dark matter haloes with masses in 



where the Cartesian components of the position vectors 
have been written (xi, z), with x^ = (x, y). Thus s z and 
u z specify the z-components of the redshift space position 
vector s and the comoving peculiar velocity field u, scaled 
in units of the Hubble parameter, respectively. Note that 
we take u to be negative for convenience. 



B. Higher order correlations 

We may now compute the correlation hierarchy for 
such a distribution of tracer objects. For a definition 
of the higher-order clustering statistics in configuration 
space and their Fourier space dual counterparts we refer 
to Appendix [X] There may also be found useful symme- 
try properties that we exploit throughout. 

Following [2(| EH, H2, |6(j, the real-space 3-point 
correlation function in the halo model, for dark mat- 
ter, haloes or galaxies, is the sum of three terms: the first 
represents the case where all three points in space are 
contained in a single halo; the second is the case where 
two points are located in one halo and the third is in a 
separate halo; the third is the case where three points are 
located in three distinct haloes - we shall refer to these 
as the 1-, 2- and 3-Halo terms and (C« )1H ) C,2H> C,3h)- 
These are written: 



(4) 
(5) 

(6) 
(7) 

I 

the range M to M + dM . 

The inverse Fourier transforms of these 3-point func- 
tions are the redshift space bispectra (c.f. Eq. IA18[) . 
They are written: 



= Ca,lH( S l: S 2, s 3) + C,2h( S 1 J S 2,S 3 ) + Cq,3h( S 1> S 2,S 3 ) J 

= i J dM d 3 y [Wy 3 n(M) J[ { U s (y - Si \M)] ; 

= 4 / II {dMidWW^MMi) U s {ji - Bi\Mi)\ [W a ]i l/'(yi-S3|Mi) 

Pa J .={1,2} 

x ^ c (yi,y 2 |M 1 ,M 2 )+cyc ; 
= ^Jf[{dM t d 3 Vl [W a ]i n(Mi) U s ( Yl - SilMi)} ^ c (y 1 ,y 2 ,y 3 \M 1 ,M 2 , M 3 ) 
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B'(ki,k aj ka) 
B' ilH (ki,k 2 ,ka) 

B'anCki.ka.ka) 



BJ )1H (ki,k a , k 3 ) + ^ 2H (kx, ka, k 3 ) + B'^OcLka, k 3 

3 



(8) 

^ y dAf[T^ Q ] 3 n(M)[]{ tf'(k,|M)} , (9) 
4 / II { dM * ^ n(Mi)C/ 5 (ki|Mi)} [W Q ]i C/ s (k 3 |Mi)P h s c (k 2 |Mi, M 2 ) + eye , (10) 



i={l,2} 
3 



B; i3H (ki,k 2 , ka) = jpJf[{dMi [W a ]i n(M t ) U'Qt&Mi)} B^fa,]^, k 3 |A/i, M 2 , 



M 3 



(11) 



where P^ c and £?g c are the Fourier transforms of the 2- 
and 3-point halo center correlation functions. 

Several advantages are gained from transforming to 
Fourier space. Firstly, once the integrals over mass are 
included, £Jh requires evaluation of a 12-D integral — the 
corresponding term _B| H is significantly simpler. Indeed, 
for the case of real space - not redshift space - and for 
spherical haloes, it is possible to write J3 3 h as the prod- 
uct of 3 2-D integrals. The calculation is slightly more 
complicated in redshift space but, as we show below, it 
remains tractable. Thus, to compute the redshift space 
power spectrum and bispectrum, we require three com- 
ponents: the abundance of dark matter haloes n(M); a 
model for the redshift space density profile; and a model 
for the inter-clustering of dark matter haloes in redshift 
space. In the following sections we describe our choices 
for these quantities. 



III. INGREDIENTS 



A. Halo abundances and bias factors 



r 



The halo bias factors associated with this mass function 
are reported in we use these in what follows. [38[ 
describe other empirical approaches to determining halo 
bias parameters. 



B. Density profiles in redshift space 

Consider the 6-D phase space density distribution 
function for dark matter particles within a particular 
halo, denoted .F(x, u|M). The density profile and veloc- 
ity distribution function may be obtained by marginaliz- 
ing over velocities and positions, respectively: 

p(x)=M J du«F(x,u|M) ; V 3D (u) = J dxF(x,u\M), 

where M is the normalizing mass. The redshift space 
density profile can be readily obtained from the phase 
space distribution through transformation to the new 
random variable s, given by our fundamental mapping 
(Eq.EJ). Hence 



The halo mass function n(M) plays a central role in 
the halo model. It has been the subject of much detailed 
study [fill, [H, [6^, • These studies suggest that, in ap- 
propriately scaled units, halo abundances should be ap- 
proximately independent of cosmology, power spectrum 
and redshift. These models for n(M) also predict that the 
real-space clustering of halos should be biased relative to 
that of the dark matter [63] ; the way in which real-space 
halo bias depends on halo mass is related to the shape of 
n(M). Thus, once the mass function has been specified, 
the problem of describing halo clustering reduces to one 
of describing the clustering of the dark matter. Of the 
many recent parametrizations of n(M), [H, HH, 16(1 loTf. 
we use that of Sheth & Tormen [63[ . Changing to that of 
Warren et al. [66| for instance, does not affect the large- 
scale matter predictions, and changes the results in the 
non-linear regime by a few percent. Note that for the bis- 
pectrum, as was shown by [521 ] , a more important issue to 
be aware of is the finite volume effect, which can signifi- 
cantly change the measured statistics for small volumes. 



p s (s±,s z ) = M J dz dx± du z du± J r (x±, z,u±,u z \M) 
xS D (s z - z + u z ) 5 D (s± - xj_), 
= M J dz du± J-(s±, z, u±, z — s z \M) . (13) 

We now assume that the density distribution of mat- 
ter within each halo is well described by a spherically 
symmetric density profile and that the particle orbits 
are isotropic and independent of position within the 
halo. Thus, the phase space distribution is separable, 
i.e. F(x,u\M) = p(x\M)Vsd(u\M)/M. Since the ve- 
locity distribution function is isotropic it may now be 
written as the product of three independent distribu- 
tions in the three coordinate directions: V 3 d(u|M) = 
Vm(ux\M)ViD(uy\M)ViD(u z \M); and we shall hereafter 
use the notation that Vid = V. Hence, 

p s (s ± ,s z \M)= [ dzp(s x ,z\M)V(s z -z\M). (14) 
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Fourier transforming yields the compact expression 

U"(k±,k z \M) = U(k\M)V(fik\M) , (15) 

where kj_ denotes a 2-D wavevector perpendicular to the 
distortion, and k z = fik = k • z is parallel to it. 

This expression shows that the redshift space profile is 
anisotropic because the spherically symmetric real-space 
profile U has been convolved along the line-of-sight di- 
rection with displacements generated by the velocity dis- 
tribution V. That is to say, V is the quantity in the 
model which generates FOG distortions and since it rep- 
resents virial motions, it is clearly the sort of nonlinear 
effect that PTs based approaches must model 'ad hoc'. 
We also note that if V makes the redshift profile substan- 
tially anisotropic, then these nonlinear effects may extend 
farther into the linear regime than one might otherwise 
have expected. 

We caution that this simple model of the halo phase 
space will almost certainly not be strictly valid, since 
it clearly neglects many aspects of the more complex 
physics that we understand to play an important role for 
the internal dynamics of haloes. [97| Nevertheless, results 
from numerical simulations show that the isotropic model 
is a good approximation, which allows one to write down 
simple expressions that provide physical insight. [98[ 

Precise details of the models we employ for the den- 
sity profile of dark matter and for the 1-Point distribu- 
tion function of velocities are presented in Appendix [B] 
Note also that owing to these models employing different 
conventions for the halo mass we must convert between 
them, and we do this using our procedure from [72| . see 
Appendix IB 21 for some details. 



C. Halo center clustering in redshift space 

On large scales the success of our analytic model will 
primarily be determined by its ability to reproduce the 
large-scale clustering of the halo centers. For this we use 
the redshift space Halo-PT developed in [48[, which is 
accurate up to the 1-loop level in perturbation theory. 
The main result we draw from that work is the idea that 
the halo density field may be written as perturbation 
series that involves the standard density PT kernels @ 
and the non-linear bias parameters Q. Explicitly we 
have the series, 



6 s hc (k,a\M,R) = £PiW] B Phc(k|M,ii)] B (i 6 ) 

n=0 

[6Um,R)] n = /n{(^3^(%)}(27r) 3 [^(k)] n 
x Z% s (q 1 ,...,q n \M,R), (17) 



of the velocity field. The functions Z„(qi, ...,q„|M, R) 
are the redshift space Halo-PT kernels symmetrized in all 
of their arguments and we make explicit their dependen- 
cies on halo mass and the scale over which the density 
field has been smoothed. Kernels up to second order are 



ryhc 

x - 

z\ c = 

ryhc 

^1.2 — 



^o hc = bo ; 

Fl lc + /(nv? Gi [l + bo] ; 

Fi% + /(fi)Mi 2 Gi, 2 
Z I q 1 L 



(18) 
(19) 



+^G 2 

92 



F-, 



> 1C + f(£l)ix\G\ 



+ i [/(Sl)/^i 2 fci 2 ] 2 — — GiG 2 b , 

where we have adopted the short-hand notation: 

= Z* c (q il ,...,q in \M,R) ; 
= Fh c (q n ,...,qjM,i?) ; 



(20) 



71 lie 



F, 



Gi 



and where 



w(|q» 



f-^ii ■ --in 



|q<i 



•QiJ-^GnCqi!,. 



(q<i H 1- q»„) • z 



(21) 
(22) 



The quantities F^ c are the nth order Halo-PT kernels 
(see Appendix [Cland |38j for complete details). The func- 
tions G n represent the nth order Eulerian PT kernels for 
the divergence of the velocity field Note that these 
expressions are almost identical to the redshift space PT 
kernels derived by [35l . |37| , however they differ in some 
subtle ways: one, we have explicitly included their de- 
pendence on the smoothing filter W(q), which is needed 
to facilitate the Taylor expansion; and two, we are ap- 
plying this in the context of haloes and not galaxies and 
so they depend on the non-linear halo bias parameters 
bi(M), instead of the nonlinear galaxy bias parameters 
(see discussion in Section IIII A|) . Note that we have also 
included bo(M), since this does not have to be zero, al- 
though we will take it to be so for all our later analysis. 

Following standard methods for calculating polyspec- 
tra, we find that the halo center bispectrum, B^ c 123 = 
B^ c (ki, k 2 , k 3 |Afi, M 2 , M 3 , R), up to fourth order in red- 
shift space Halo-PT, is 

£W23 = 2P 11 (k 1 )P 11 (k 2 )Z i l c (k 1 \M 1 ,R)Z^(k 2 \M 2 ,R) 
x Z 2 hc (k 1 ,k 2 |M 3 ,i?) + 2cyc . (23) 



where D(t) is the linear theory growth function and 
/(f2) = d log D{a)/d log a is the logarithmic growth rate 



Inserting our expressions for the redshift space Halo-PT 
kernels, Eqs (fT8l 120]) . in to the above expression, reveals: 
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W(k 1 R)W(k 2 R)W(k 3 R) 

x F 2 (k 1 ,k 2 )+/3 3 /i?2G 2 (ki,k 2 



= 2b 1 (M 1 )b 1 (M 2 )b 1 (M 3 )P 11 (k 1 )P 11 (k 2 ) 11 {l + pnj$[l + bo(Mi)]} 



W(hR)W(k 2 R) fc 2 (M 3 ) 1 a Ml ^ 2 



5/(n)Mi2fcu ^ [1 + /3 3 A4] + g [1 + & Mi] 



M2 



2 eye. 



(24) 



where fa = /(f2)/&i(Mj) and where Cj(Mi) — 
bj{Mi)/b\{Mi). As in our real space work on the power 
spectrum, we are now faced with the situation that we 
have solved for the bispectrum of halo centers filtered on 
scale R, and in fact we would like to recover the unfiltered 
bispectrum. As in [38j , we take the following ansatz: the 
filtering of the spectra can be reversed through the fol- 
lowing operation: 



P(k) 
£(ki,k 2 ,k 3 ) 



P(MR) . 

W 2 {kR) ' 

B(ki,k 2 ,k3|fl) 

W(k 1 R)W(k 2 R)W(k 3 R) 



(25) 



and this explains the form of the left-hand-side of 
Eq. ([24]) . An alternative approach to the filtering issue 
for the power spectrum and bispectrum was proposed by 
[73 | , we shall leave the solution of this problem for future 
consideration. 

With these ingredients prepared, we are now in full 
possession of a complete description of the bispectrum of 
galaxies, haloes and dark matter in the halo model and in 
the presence of a local, non-linear scale dependent bias. 
In the next section we develop these equations further. 



IV. THE BISPECTRUM MONOPOLE 

A. Euler angle averages 

The redshift space bispectrum is an anisotropic func- 
tion on the sphere that depends on 5 variables. The first 
three are the triangle configuration, and the other two 
specify its orientation with respect to the z-axis. How- 
ever it is common practice to measure this quantity av- 
eraged over all possible orientations of the coordinate 
frame. Thus to compare with this isotropized observ- 
able, we shall now derive the isotropized form for the halo 
model, which we shall denote B . Interestingly, the fol- 
lowing approach is identical to that which one performs 
for the triaxial halo model in real space [20|, |72| , since on 
small scales, one may effectively think of transforming to 
redshift space as simply transforming a set of spherical 
haloes into a set of prolate ellipsoids whose semi-major 
axes all point along the line-of-sight. Of course, the real 
situation is more complex, since the halo centers are also 



distorted according to the halo velocity projected along 
the line-of-sight, but nevertheless we may borrow some 
mathematical machinery from the triaxial halo analysis. 
The isotropic function is thus, 



B (k u k 2 ,9 12 ) 



1 



CZ71 d(cos7 2 ) d73 



x B s [H(m, 72, 73)ki,ft(7i, 72, 73)k 2 ] ,(26) 

where 72. (71, 72, 73) is the rotation matrix for the compo- 
nents of the basis vectors. 72.(71,72,73) is parametrized 
by three position or Euler angles, (71,72,73) and these 
are the z' — y' — z" rotation angles (see Appendix [D] for 
the explicit form of the matrix we use). Note that we 
assign uniform probability on the sphere to each triple of 
angles. 

The following considerations simplify this expression 
considerably. Firstly, in the PT expressions for B s , each 
term depends on either the angle between two k-vectors 
or the projection of each vector along the z-axis. In the 
first case, the use of matrix notation shows that 

pcjfkj = [72(7i, 72, 7 3 )k»] T 72(7i, 72, 7 3 )k, 
= kf72(7i,72,73) T 72(7i,72,73)kj 
= kfk; ; (27) 

this is the well-known result that scalar products are in- 
variant under rotations of the coordinate basis functions. 
In the second case, the projection of each rotated vector 
onto the line of sight direction (or z-axis) can be written 

k' • z = 72(7i, 72, 73)k • z 

= k x sin 72 cos 71 + k y sin 7 2 sin 71 + k z cos 7 2 

= A(k, 7 i,72) , (28) 

Thus, we see that the resultant function must be invari- 
ant under the 73 rotation, and hence this integral may 
be computed trivially. Finally, because our final quan- 
tity B must be independent of the initial locations of the 
fc-vector triple, we may without loss of generality choose 
these locations to be as convenient as possible. There- 
fore, we let the initial ki vector lie along the polar-axis 
and constrain k 2 and k 3 to lie in the z—y plane: i.e. 

[ki] T = (0,0, kx) 

[k 2 ] T = (0, fc 2 cos#i2, fc2sin0i 2 ) 

[k 3 ] T = (0, -fc 2 cos6»i 2 , -fci - fc 2 sin(9i 2 ) (29) 
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where the last equality uses the closure condition: 
k 4 = 0. Thus, the fc-vectors rotated into the new 
basis and dotted with the z-direction are now written: 

A a = cos 72 (ki) z ; 

A 2 = SUI72 sin 71 (k 2 ) y + COS72 (k 2 ) z 5 

A 3 = -sin 72 sin 71 (fe 2 )„ - COS72 [{ki) z + {k 2 ) z } .(30) 

The parameters A4 = A(k i ,7 1 ,7 2 ) are simply related to 
the cosines of the fc-vectors along the z-axis: 

Ax A 2 kx k 2 

Mi = 7- ; M2 = 7- ; fj- 3 = -fj-i-, — M27-, (31) 

kx k 2 k 3 k 3 
I 



where q = k 2 /kx- 



We may now apply the above operation directly to our 
expressions for the anisotropic bispectrum (Eqs [9l- lll[) . 
On inserting the A4 into each instance of fa in the density 
profiles (Eq. [T5|) and the halo center power spectra and 
bispectrum (Eqs IB"4l k, [24]). we find that the 1-, 2- and 
3-Halo terms for the bispectrum monopole become 



B m (k u k 2 ,9 12 ) 



d 7l d(cos7 2 ) dM[W a ] 3 n(M)l[{U{ki\M)V(faki\M)} ; 



(32) 



B S 2l i{kx, k 2 ,0x 2 ) 



={1,2} 



d 7l d(cos7 2 ) J] { dMilW^ibxiMAniMAUikilMAVifiikilMi) 



:[W a ]x U(k 3 \Mx) V{^h\Mx)Pxx{k 2 ) { 1 + [A + ft] £ + ft ft M || + cyc . (33) 



B 3H (fci, A; 2 ,^i2) 



1 



i=l 



d7id(cos 7 2) Hi / dMilWaliniMAUikilM^VifakilM^ 



x B£ c (fcx , k 2 , e 12 , Mi , M2 1 Mi , M 2 , Af 3 ) . 

I 



(34) 



The only variables that depend on the Euler angles 71 
and 72 are fa and A;. Each term requires evaluation of 
no more than a 4-D embedded integral: two integrals for 
the Euler angles, one for the mass, and one for the Fourier 
transform of the density profile. Note that for simplicity, 
we have kept only the leading order contribution to the 
the 2-Halo term. Technically this should be taken up 
to the 1-Loop level to be consistent with the bispectrum 
which is 4th order in 5. However, this issue is beyond the 
scope of the current paper and will be addressed in 



B. Computational considerations 



Our expressions for the bispectrum as presented above 
are complete. However some calculational effort is still 
required before we may attempt a practical implementa- 
tion on the computer. We now present some simplifica- 
tions. 

We begin by defining some convenient notation: Let 



(i) 

ip s j and ip v j denote the following integrals: 



ki,...,k,-) = 4 fdMn(M)b i (M) 

Pa J 

j 

x i[{[W a ]U{k l \M)V{^ki\M)j ; 



i=i 

Vv,i(ki, . . .,kj) 
j 



dMn(M) 



(35) 



f[{ \W a ]U(ki\M)V(inki\M)} 
1=1 



(36) 



The first integral generalizes the halo bias weighting 
scheme applied to the density field for the situation where 
j-points are within a single halo. The second integral gen- 
eralizes the weighting scheme to the similar situation for 
the halo velocity field. This notation has some similari- 
ties with that of [46[ , but is different in the way in which 
the velocity field is treated - recall that the velocity field 
has been assumed to be unbiased. 

In this notation we may re-write the 2- and 3-Halo 
terms in the bispectrum; the 1-Halo term requires no 
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simplification. Through rearrangement of the mass inte- through substitution of the appropriate kernels, we find 
grals and expansion of the halo center power spectrum that the 2-Halo term can now be written 



^, 2 (k;,k' 3 )V>g(k 2 ) + V^ +2cyc} . (37) 



The 3-Halo term is significantly more complex, owing to we may again isolate the integrals over mass and use the 
the halo center bispectrum being the product of two first ip functions to obtain 
order kernels and one second order kernel. Nevertheless, 



Ssh^i^M = ^ f d 7 id(cos 72 ) {2F(fc 1 )P(fc 2 )T 1 (k' 1 ,k 2 ) 



T 2 (k[,k 2 \k' 3 )4%' 3 ) + T 3 (k' 1 ,k 2 \k' 3 )ip v , 1 (k 3 ) + Wi2, : 



where we have defined the following useful quantities: 



T4(k' 1 ,k 2 |k 3 )^° ) (k'3 



T 2 (ki,k,-|ki) 
T 3 (k,,k,|kO 
T 4 (ki,kj|k J ) 



n {4% m )+vii[^Ak m )+m)4%rn)]} ; 



F 2 (ki,k.,) + -WijjfityfHjkij 



£ t\i<i fx, j 



Mi 2 i ft 2 



+ 2 eye | ;(38) 

(39) 

(40) 

(41) 
(42) 



r 



and 



W(hR)W(kjR) 



(43) 



W(hR) 

The advantage of this reformulation of the 2- and 3- 
Halo terms is that we have decomposed the integrand into 
a set of algebraic functions of the ip integrals (Eqs [35l - l3"6|) . 
and auxiliary functions, and these may all be computed 
in parallel making the computation highly modular. 



C. The large-scale limit 

Next we consider the redshift space bispectrum in the 
very large scale limit, as this should asymptotically re- 
duce to the standard PT expressions for the bispectrum, 
modulo discreteness corrections for the point process as- 
sociated with the halo field. However, let us first exam- 
ine the large-scale limit of Eqs (|35|) and (f3T>!) . On letting 



{ki} — > 0, the density profile terms become U s (ki) — > 1, 
the window function becomes Wijj — * 1, and so 

f = (bj(M)[W a ] j ) 

{fe.,}-o^ ([ Wa }Y 



lim ib v j 



m)([w a ] j ) 

([w a ]) j 



(44) 



where (...) = / dMp(M) with p(M) = n(M)/n H , 
njj being the total number density of haloes in the re- 
quired mass range. When j = 1 we write these functions 
more simply as: the average non-linear bias parameter 

(i) — 

for the tracer particles, ip s [ = b a ,i, and the logarithmic 
growth factor for the linear velocity field ip Vt \ = / (0). 

The bispectrum in the large-scale limit can now be 
computed directly. The 1-Halo term is trivially obtained, 
and the 2- and 3-Halo terms can be developed through 
replacing the ip functions in the general expressions (|37|) 
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and ([38)) , for their large scale forms (|44|) . After some 
algebraic manipulation we arrive at the result 



B a (k 1} k 2 ,9 12 ) = B aPT (k 1 ,k 2 ,6 12 ) 

[P(k 2 ) + P(h) + P(fci)] + 



;(45) 



1 Q,1H 



where the large scale PT bispectrum monopole in redshift 
space is given by 



Bl PT (k u k 2 ,9 12 ) = i- f d 7l d(cos 72 ) J] {b aA +fjlf(n) [l + b afi ]} 

A — 11 



x |2P(fc 1 )P(fc 2 



6 aa T 2 (k 1 ,k 2 |k 3 ) + /(n)T 3 (ki s k a |k3) + ^ + 6 Q ,oT4(ki,k 2 |k 3 ; 



2 eye ; (46) 



r 



this is equivalent to that found by [3a, |37|. We also 
defined the 1- and 2-Halo 'effective ' number densities to 
be: 



1 



"a,lH 



(w a ) 3 



(47) 



we show some tentative evidence for these effects in the 
measurements from our numerical simulations. 

Before continuing, it should also be noted that on set- 
ting /(f2) = 0, one may recover the f-Loop PT bispec- 
trum in real space for a set of biased tracer particles a. 



(bi(M)[W a } 2 ) - /(f])2 (\W a f) 



/(«) 



(w a y 

{[W a ] 2 ) 



-b a 



5 {W, 
(6i(Af)[W- a ] 2 )" 
(W a ) 3 



(48) 



Our final expression for the large-scale limit (Eq. [43]) 
is similar to the standard theoretical expectation for the 
bispectrum recovered from a Poisson point process sam- 
pling of a continuous field (c.f. Section 43 in [76j). How- 
ever, the halo model effective shot-noise terms (the last 
two terms on the right hand side of Eq. |4"5]) . are very 
different from the standard form, which would simply 
have e.g. ri aj2 H = n a . These effective number densities, 
Eqs (|4T[) and (|48|) . represent the fact that in the halo 
model we assume that dark matter haloes are Poisson 
sampled into the density field and that the tracer parti- 
cles are injected into these haloes, and so simply act as 
different weights. The issue of sampling tracer particles 
into the density field has some important implications for 
how one should extract information from galaxy surveys. 
We shall reserve this investigation for future work. 

In Appendix [F] we present a short discussion of how 
these shot noise corrections can impact the reduced bis- 
pectrum. The main results are as follows: For stan- 
dard shot noise, and in the low-sampling limit there is 
no configuration dependence and Q d = 1/3, sub-script d 
denotes discrete. For the case where shot noise is sub- 
dominant, the Q d is reduced relative to the continuum 
limit Q for all configurations. In the halo model, and in 
the low sampling limit, again there is no configuration 

dependence and Q HM = ([W a ] 3 ) {[W a ]) /3 ([W a ] 2 f '. For 
the case of sub-dominant shot noise, Q HM is not neces- 
sarily smaller than the continuum limit case. In Sec. IVII 



D. The small-scale limit and hierarchical models 

On small scales, the bispectrum is dominated by the 
1-Halo term, given by Eq. (|3"2"|) . Our understanding of 
its behavior in this limit can be guided by considering 
the case where all haloes are of the same mass. In this 
situation we have n(M) — ► 5 D (M — M')nn and p a — > 
[Wa]. Applying these conditions to Eq. ([3"2"]> . we have: 



B m (k u k 2 ,0 12 \M) 



1 



d^fi d(cos7 2 ) 



Y[{U(ki\M)VOMki\M)} • (49) 



We may follow this same procedure for the power spec- 
trum (see Eq. IE2[) and so construct the reduced bispec- 
trum, whence 



Q s 



j d 7l d(cos 72 ) nti {vj^hm/ujhiM)}/^ 

nfl( ai )n^l(a 2 )/U(k 3 \M) 2 + 2cyc 

(50) 

where TZi 2 {a\) is given by Eq. (|E9[) . Notice that this ex- 
pression no longer depends on the weights [W a ] or num- 
ber densities of tracers nn, but simply the real space 
profile and the 1-PT velocity profile. If the density and 
velocity profiles were mass independent, then Eq. (|50[) 
would predict the same configuration dependence for all 
haloes. However, for realistic redshift space profiles, this 
is not the case (see Appendix [B| . Thus, if the halo model 
is a good description for small scale clustering, then the 
hierarchical model is unlikely to be correct in real or red- 
shift space, and we may generally extend this statement 
to any tracers of the density field. Finally, and somewhat 
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interestingly, notice that if the ratio V(fiiki\M)/U(ki\M) 
is mass independent, then the configuration dependence 
of the bispectrum becomes universal, modulo an ampli- 
tude off-set. 



E. The White-Seljak Approximation 

The White-Seljak approximation (hereafter WS) is the 
supposition that haloes are randomly oriented relative to 
each other so that in computing the monopole of the bis- 
pectrum, orientation averages can be taken separately 
over each individual halo and over the large scale orien- 
tation of the halo relative to each other jia, |46| . In the 
triaxial halo model, showed that the overall contribu- 
tion from halo alignment to the matter power spectrum 
was negligible. Thus we may similarly assume that on 
large-scales the 2- and 3-Halo terms will not be sensitive 
to the orientation of the FOGs - and this allows us to 
use the isotropic redshift space density profiles instead of 
the anisotropic profiles. However, since the bispectrum is 
more sensitive than the power spectrum to the shapes of 
structure, we shall be a little more cautious, and demon- 
strate the validity of this approximation in Section IV Al 

In the WS approximation we therefore take, 



1 

47T 



d7id(cos7 2 )V'o( k 'i) 



^ 2 (k u k 2 ,9 12 ) = ^ J d 7l d(co S 72)Vg(ki,k' 2 ) ,(51) 



and with similar expressions for the tp v j functions. On 
replacement of these terms into Eq. (|37j) we find that the 
2-Halo term simplifies to: 

B ajH( fc l: fc 2,6»12) = P(k 2 ) { ( fc l > fc 3 , 013)^1 ( fc 2 ) 

+ ^ ^,2(^1^3,^13)^1(^2) 
+ ^2(^1,^3,^13)^,1(^2) 

+ r$ o ,3(to>*S>0i3)?' o ,i(*a) } +2 eye .(52) 



Similarly, on applying the WS approximation to Eq. 
the 3-Halo term reduces to: 

#S S (fci,fe,M = 2P{k l )P{k 2 ) J ^rf(cos 72 ) 



xT^kf,) T 2 (k[y 2 \k' 3 )^hk 3 ) 



+T 3 (ki,k' 2 |k^„ il (A; 3 ) + Wi2, 3 



T4(k' 1) k 2 |k 3 )^° ) (fc 3 ) +2cyc; (53) 



where 

Tl(ki,kj-) = JJ \jP S \( k m) +Mm ^,l( fc m) 

+ f(n)^f}(k m )]} . (54) 

This completes our analytic investigation. In the next 
sections we shall provide numerical evaluation of our ex- 
pressions for the bispectrum. 

V. ANALYTIC RESULTS 
A. Testing the WS approximation 

Figure [1] compares the predictions for the reduced 
redshift space matter bispectrum (Eqs (f32|) . ([37]) . and 
([38]) ) with the WS approximate expressions Eqs (j32j) . 
(IS"2")) . (T55|) . The four panels show fc-space triangles with 
k 2 /ki = 2 and with fci = {0.05, 0.1, 0.5, 1.0} ft Mpc -1 . 

On very large scales (top left panel), k = 0.05ft. Mpc - , 
the predictions are dominated by the 3-Halo term, and 
the approximate expressions (thin blue lines) are in al- 
most perfect agreement with the exact Halo Model pre- 
dictions (thick red lines), with < 2% deviations across 
different configurations. This is to be expected follow- 
ing our derivation of the bispectrum in the large-scale 
limit, c.f. Eq. (f4"5"|) ; i.e., there is no dependence on the 
halo profiles and so the WS approximation does not play 
a role here. On slightly smaller scales (top right panel) 
k\ = 0.1ft. Mpc -1 , small but significant departures are 
found at the level of < 10%. It can be seen that these 
are entirely due to deviations in the 3-Halo term - the 
order in which we spherically average the profiles is now 
important. However, on still smaller scales (bottom left 
panel) fci = 0.1ft Mpc -1 , whilst the discrepancy between 
the approximate and exact 3-Halo term becomes larger, 
the 1- and 2-Halo terms begin to dominate and so the 
difference in the total appears < 6%. Finally, on the 
smallest scales considered, the 1-Halo term comes to fully 
dominate and since no approximation is made here the 
results are in good agreement, < 5%. 

For a given scale, the largest deviations typically ap- 
pear for the case of colinear triangles, i.e. where all 
three fc-vectors are co-linear. This leads us to suppose 
that the equilateral bispectrum as a function of scale 
will show agreement at the level of < 5% across a wide 
range of scales under this approximation. Since current 
observational measurements of the bispectrum on large 
scales have sample and cosmic variance errors of the order 
~ 50% on scales k\ ~ 0.1ft Mpc -1 , going down to several 
percent at k\ ~ 1ft Mpc -1 , we anticipate that, at least 
for current data, the WS approximation should be use- 
ful. We highlight again that the main advantage of this 
approximation is the increased speed with which one can 
compute the bispectrum: the tpg^- and ip v j functions are 
only evaluated once for a particular configuration, as op- 
posed to thousands of times. However, in all that follows 
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FIG. 1: Configuration dependence of the reduced bispectrum monopole in redshift space - comparison of the predictions from 
the exact Halo Model expressions with those under the WS-approximation. Top sections of each panel show results for triangles 
with k2/ki = 2 and on scales: ki = {0.05, 0.1, 0.5, 1.0} h Mpc" 1 . The dash, dot-dash and dotted lines in each panel correspond 
to the 1-, 2- and 3-Halo terms, respectively. The solid lines correspond to the sum. Thick (red) and thin (blue) lines are the 
exact Halo Model and WS approximate results, respectively. Bottom sections of each panel show the ratio of the total WS 
approximate bispectra to the exact Halo Model calculation. Line styles have same meaning as in top panels. 



we shall only show results for the exact evaluation (no 
WS approximation) of our redshift space Halo Model. 

B. Comparison with perturbation theory 

Figure fj] compares the analytic predictions for the re- 
duced bispectrum from our model with corresponding re- 
sults from PT. The four panels show again A:-space trian- 
gles with with k^/ki — 2 and for the same scales as pre- 
sented in Fig. [TJ The solid (red) lines in each panel show 
our HMs predictions (recall HMs means Halo Model in 
redshift space) in the WS approximation. The (red) dash 
lines show the PTs predictions, as given by our Eq. (|46[) . 
The (blue) dotted lines correspond to real space PT pre- 



dictions. 

On the largest scales k\ = 0.05 /iMpc -1 (top-left 
panel), the HMs and PT results match almost perfectly: 
the configuration dependence shows excess signal for col- 
inear triangles, indicating that on these very large scales 
non-linearity induces structures that are, on average, 
more filamentary than spherical [35|, [75[ (in this diagram 
spherical perturbations are best probed by isosceles tri- 
angles, 6*i2 ~ 27r/3). However, we notice that there is 
a small deviation < 5% for the situation where the k\ 
and &2 vectors are aligned. This owes to the fact that 
the 1- and 2-Halo terms are non- vanishing as k — > 0, 
and as discussed in Section IIV C[ this gives rise to an 
'effective shot-noise' like behavior. A discussion of why 
co-linear k\-ki triangles (i.e. 0\2 — 0) are more prefer- 
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FIG. 2: Comparison of analytic model predictions with results from PT. The configuration dependence is shown for for triangles 
with fe/fci = 2 and on scales: k\ = {0.05, 0.1, 0.5, 1.0}/i Mpc -1 . In all panels: solid (red) lines represent redshift space halo 
model predictions (HMs); dash (red) lines correspond to redshift space tree- level PT (PTs) [35j; and dotted (blue) lines 
correspond to real space tree-level PT. 



entially affected is given in Appendix [F]) . It should also 
be noted that unlike for the case of standard shot noise, 
Q HM > Q FT for configurations close to isosceles trian- 
gles. 

On slightly smaller scales k\ ~ 0.1 ft Mpc -1 , both the 
PT and PTs predictions show a small increase in config- 
uration dependence, with the PTs having slightly more 
signal for colinear triangles than PT. The HMs predic- 
tions are in qualitative agreement with PT. However, the 
flattening off seen in the previous panel is now much more 
apparent. Recalling the corresponding panel in Fig. [1] it 
can be seen that this is attributed to the rapidly rising 
1- and 2-Halo terms. 

On smaller scales still k\ = 0.5 ft, Mpc -1 , the PT and 
PTs predictions continue their previous trends, exhibit- 
ing a slightly increased configuration dependence. How- 
ever, the HMs predictions are quite different, having a 
very strong [/-shaped configuration dependence. This 
owes to the 1- and 2-Halo terms becoming dominant. 

On the smallest scales considered fci = 1.0 ft Mpc -1 , 
the HMs predictions show a dramatic configuration de- 



pendence, with a very strong signal for colinear triangles 
and a broad plateau for triangle shapes around isosce- 
les configurations - this is the ' [/-shape', which was first 
noted in the bispectrum by [35j and in the 3-point cor- 
relation function by [36[ ■ Recalling Fig. Q] (bottom left 
panel) , we see here that the prediction is completely dom- 
inated by the 1-Halo term, and so this [/-shape feature 
is simply an imprint of the halo shape in redshift space 
- the FOG. This result constitutes a more direct demon- 
stration of the discussion from Section IIVD| that hier- 
archical models are unlikely to be a good description for 
the higher order clustering statistics. 

In this section we have shown from purely theoreti- 
cal considerations, that to use the galaxy bispectrum on 
large scales k\ < O.lftMpc -1 as a precise tool for cos- 
mology, one needs to understand exactly how to include 
the FOG effect into the modeling and also how to include 
non-trivial discreteness effects of matter. In the next sec- 
tion we confront the model with results from numerical 
simulations. 
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VI. COMPARISON WITH NUMERICAL 
SIMULATIONS 

A. Numerical simulations 

In order to test our redshift space bispectrum we gen- 
erated an ensemble of 8 LCDM simulations, these were 
identical in every way, except that for each simulation 
different random realizations of the initial Fourier modes 
were drawn. The cosmological parameters for the en- 
semble were selected to be in broad agreement with the 
WMAP best fit model [|: O m = 0.27, Q A = 0.73, 
£l b = 0.04MI = 0.72 and a 8 (z = 0) = 0.9. We used the 
cmbfast [571 ] code to generate the linear theory trans- 
fer function, and we adopted the standard parameter 
choices, and took the transfer function output redshift 
to be at z = 0. The initial conditions for each simula- 
tion were then laid down at z = 49 usin g the p ublicly 
available 2LPT initial conditions generator [73, [78j] . Sub- 
sequent gravitational evolution of the equations of motion 
was then performed using the publicly available Gadget 2 
code |79j . Each simulation was run with N — 400 3 par- 
ticles in a comoving volume of length L = 512 3 Mpcft. _1 
and with a comoving force softening set to 70 kpc//i. 
The simulations were performed using a 4-way dual core 
Opteron processor system, and each ran to completion in 
roughly ^1800 timesteps from redshift z = 49 to z = 
and this translates to ~2 days of wall clock time. 



B. Estimating the power spectrum 

Before comparing our analytic model with the bispec- 
tra estimates from the simulations, it is instructive to 
compute the real and redshift space power spectra of the 
z = outputs of the ensemble. 

The density Fourier modes are estimated using the 
conventional Fast Fourier Transform (FFT) method: the 
dark matter particles were assigned to a regular cubical 
grid using the 'Cloud-In-Cell' (CIC) scheme gfjj. The 
FFT of the gridded density field was then computed us- 
ing the publicly available FFTW routines [8l| . Each result- 
ing Fourier mode was then corrected for the convolution 
with the mesh by dividing out the Fourier transform of 
the mass-assignment window function. For the CIC al- 
gorithm this corresponds to the following operation: 



<5 d (k) = <5 g (k)/W C ic(k) 



where 



Wcic(k) = n 



2 = 1,3 



sin [7rfci/2/cNy] 

[7Tfc l /2fcNy] 



(55) 



(56) 



where sub-script d and g denote discrete and grid quan- 
tities, and where k-^y = nN g / L is the Nyquist frequency 
of the mesh and N„ is the number of grid cells. 




k [h Mpc : ] 

FIG. 3: Real and redshift space power spectrum of dark mat- 
ter particles at z — measured from the ensemble of LCDM 
simulations. Top panel: The power spectrum. Blue and 
red points show real and redshift space quantities, respec- 
tively. The solid, dot-dash and dotted lines show the total 
halo model, 1-Halo and 2-Halo terms. The triple dot-dash 
curve shows the predictions for the Poisson shot noise error. 
The dash line shows the linear theory. Bottom panel: The 
measured power spectra ratioed with a no-baryon linear the- 
ory power spectrum that has the same overall transfer func- 
tion shape. Note that for the redshift space spectra, we have 
scaled out the Kaiser boost. 



The power spectra of the discrete particles on scale fc; 
are then estimated by performing the following sums, 



V, 



M 

-Y 



|<*d(kz 



(57) 



where M is the number of Fourier modes in a spherical 
shell in fc-space of thickness Afc. Note that the mode-by- 
mode correction differs from the analysis of [54],[82] where 
the correction for charge assignment was performed by 
computing the spherically averaged window and dividing 
it out. Mode- by-mode correction for the power spectrum 
was also performed in [38l |75| . 

Figure [3] shows the mean and l-cr errors for the power 
spectra of dark matter particles in both real and redshift 
space. The errors are computed directly from the 8 real- 
izations. The spectra were computed using a 1024 3 FFT 
and we show all frequencies up to the Nyquist frequency 
- the highest k- modes show signs of increased power from 
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both Poison shot noise (triple dot-dash green line) and 
the aliasing of power from smaller scales. We note that 
on the largest scales probed, the power spectra show a 
sequence of wiggles, these are the well known Baryonic 
Acoustic Oscillations (BAO) that have been discussed 
much in recent times [3, [H, [H, [H, [85[ - we shall not 
discuss these in this paper. The solid lines show the 
total halo model predictions in real (blue lines) and red- 
shift space (red lines). We see that, whilst the real space 
model does reasonably well, with an accuracy of the or- 
der ~ 10%, the redshift space predictions fare less well, 
especially for scales k > 0.2/iMpc . Here the model sys- 
tematically underpredicts the data by roughly ~ 20%. 

These predictions were very sensitive to how we mod- 
eled the FOG effects, i.e. the ID velocity dispersion of 
particles within haloes v(k\M). Originally, we had sim- 
ply used Eq. (|B12[) with e = 1, however in this case 
the predictions were particularly poor. We therefore 
investigated M — am relation in our simulations more 
closely. Using a standard Friends-of-Friends (FoF) algo- 
rithm with link length / = 0.2, we located all haloes with 
M > 1.0 x X0 13 M e /h. The ID Velocity dispersions were 
then estimated for each halo and binned as a function of 
mass. We found that e = 1 overpredicted the measured 
values by > 20%. Fitting on e it was found that e = 0.76 
provided a better fit to the data, but it was not perfect, 
it having an accuracy no better than ~ 10% (see Fig. 
As can be seen from the figure using e = 0.76 does not 
generate the correct power spectrum. 

A possible reason for this discrepancy could be that 
the 2-Halo term in the redshift space power spectrum 
only includes linear halo motions - non-linear terms do 
contribute to this term (see [H, 0, [85| ) and inclusion of 
these non-linear corrections may help alleviate this prob- 
lem. Another, is that the concentration-mass relation, 
which is vital for getting the correct normalization of the 
halo density profiles, may not be sufficiently accurate. In- 
deed we note a small discrepancy between the real space 
measurements and Halo Model predictions. Recent im- 
provements on this relation by [86j may help to alleviate 
this problem. Also, we have neglected the scatter in the 
concentration parameter, and it is well known that this 
can change the predictions by a few tens of percent in 



with some changes. Our estimator can be written: 



the nonlinear regime 87 1. Additionally, there is the is- 
sue of halo triaxiality 20, [zl] . We shall not pursue these 
subtle corrections here, since our purpose is to simply 
present the theoretical framework and show that it gives 
reasonable agreement with the simulation data. 



C. Estimating the bispectrum 

Our estimator for the bin and spherical averaged bis- 
pectrum was developed following the work of [75] , but 
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A practical implementation of this estimator involves 
computing the following sum 
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where is an integer vector from the /c-space origin to 
a mesh point and so labels the modes, and where N tT i 
represents the number of independent momentum con- 
serving k- vector triangles in the shells V\ and V 2 . In the 
above we take only the real part of the product of the 
three Fourier modes, owing to the reality of the bispec- 
trum (see Eq. IA21j) . Note that when computing the sums 
over fc-space triangles we randomly sample modes from 
the set of all possible triangles. Typically we limit the 
computations to 10 4 modes per shell (i.e. 10 s triangles). 
This method gives a sufficient number of independent k- 
space triangles for a high accuracy estimate, whilst keep- 
ing code execution times tolerable. 

To accurately estimate Q we are also required to esti- 
mate the combination Q fac = P(ki)P(k 2 )+P(k 2 )P(k3) + 
P(k3)P(ki). There are several approaches to achieving 
this: one, we could simply estimate the power spectrum 
as in Eq. (|57|) and then construct Qfac from this; alter- 
natively one can compute an estimate of Q iac using only 
those same modes that are used to estimate B. We adopt 
this latter approach since expect that it will reduce sam- 
ple variance. Along with our estimates for B, we also 
therefore accumulate 
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(62) 



Q faC:d (# 12 ) = P hd P 2 ,d + P2,dP3,a + P 3 ,aPi,d ■ (63) 

We draw close attention totiie factjjhat these estimates 
for the power spectra Pi, P 2 and P3 are not the same 
as in Eq. ([57]) : in the above case we average over all 
/c-space triangles that are used and not just the unique 
modes. We have also made it explicitly clear that the 

estimates for Pi and P 2 do not change with 8i 2 , but 

that P3 does. Following this procedure helps to re- 
duce cosmic variance. We also note the following pit- 
fall: had we estimated <5f ac for each fc-space triangle 



15 



and then aver aged these es timates over all triangles, i.e. 
taken Q fac = P(ki)P(k 2 ) + 2cyc, then we would have 
been angle averaging products of power spectra. In real 
space, where P is an isotropic function on the sphere, this 
makes no difference, however in rcdshift space, where P 
is anisotropic, this approach would be incorrect and no 
U-shape would be seen in Q. 

In order to correct B and P for discreteness we also 
estimate the shot noise terms as [761 ]: 
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(64) 
(65) 
(66) 



Estimates of the shot noise corrected continuous spectra 
are then arrived at through the following set of opera- 
tions: 



P = P d - P shot ; 

B = Bd — £>shot ; 

Qfac = Qfac, d — Qfac, shot i 

3 = t/Sfac • 



(67) 
(68) 
(69) 
(70) 



In our measurements we shall show both spectra with 
and without the shot noise corrections. 

Lastly, as a consistency check, we also estimate the 
imaginary bispectrum, which is given by Eq. (|60[) only 
now we take the imaginary piece of the product. Thus, 
if our estimate is correct, then this quantity should on 
average be zero. However, as the number of independent 
triangles becomes small the imaginary piece may become 
non-zero due to statistical fluctuations. 



D. Bispectrum results 

Figure 0] shows the mean and the l-cr errors on 
the mean (i.e. we divide the standard error bars by 
l/v7) for the configuration dependence of the bispec- 
trum B(ki, k-z, #12) in both real and redshift space for 
the 8 LCDM simulations. The four panels show re- 
sults for fc-space triangles with ki/k\ = 2 and k\ — 
{0.05, 0.1, 0.5, 1.0}/iMpc _1 and thin (blue) and thick 
(red) lines distinguish between real and redshift space 
quantities. 

Considering the largest scales (k = 0.05 /iMpc -1 ), we 
see that the ensemble estimate is rather noisy, this owes 
to the large sample variance on these scales. Neverthe- 
less, it can be discerned that the redshift space (large 
solid stars) estimate has a slightly higher amplitude than 
the real space (solid points) estimate. We also note that 
there are a few bins that give a significant non-zero con- 
tribution to the imaginary bispectra. We attribute this to 
the large sample variance errors, but in the main these 



points are below the true signal. In all cases the re- 
sults for the bispectra appear to be consistent with the 
PT and HM predictions, to within the errors. To say 
anything more definite on these scales will require much 
larger simulation volumes, and we shall defer this ques- 
tion for future study. 

On intermediate scales (fc = 0.1 /iMpc -1 ), the esti- 
mates for the real and redshift space spectra are much 
more significant, the errors being well defined. This is 
supported by the fact that the imaginary spectra have 
amplitudes that are now too small to be plotted. It 
can be seen that the results have a strong dependence 
as a function of triangle configuration. On comparing 
with the PT and PTs, we see that on these scales PT 
under-predicts the measured quantity by roughly ~ 20%, 
whereas PTs overpredicts the amplitude by a similar 
amount. Considering the predictions of the HM and 
HMs, it is clear that in real space the model is a sig- 
nificant improvement over PT; whereas in redshift space, 
whilst the configuration dependence appears to have been 
slightly improved, the amplitude is still too large. 

On smaller scales (fc = 0.5 hMpc^ 1 ), we see that the es- 
timates of the spectra are of even higher signal-to-noise 
and that they form tight-loci across the configuration. 
Again the imaginary bispectra are insignificant. It may 
be noticed that there is a large difference between the 
real space PT and the simulation estimate - almost two 
orders of magnitude, this is due to the fact that B ~ QP 2 
and the power spectrum is significantly larger than lin- 
ear. However, this dramatic change in the measured bis- 
pectrum is captured relatively well by the HM, which 
underpredicts the result by only ~ 20%. Turning to the 
redshift space results, we are surprised to see that the 
PTs is of the same amplitude and has somewhat similar 
configuration dependence as the estimate. This agree- 
ment seems coincidental, given the poor agreement in 
real space. In reality the true dynamics on these scales 
must already be very non-linear. Lastly, we draw at- 
tention to the fact that our HMs result is in excellent 
agreement with the simulation data. 

Examining scales on the order of the virial radii for 
clusters (fc = 1.0/iMpc -1 ), we see again that the es- 
timates are of very high significance. Again, the HMs 
predictions are in excellent agreement. This essentially 
vindicates our form for the 1-Halo term, and means that 
the configuration dependence is very much governed by 
the FOG distortions. 



E. Reduced bispectrum results 

Figure [5] shows similar results as in Fig. [4] but for 
the configuration dependence of the reduced bispectrum, 
Q(ki, k%, 812). Again errors are shown on the mean of Q. 

On the largest scales probed, (k = 0.05 ft-Mpc -1 ), the 
estimates are noisy, but there is evidence for an excess 
of signal for co-linear triangles - meaning that on aver- 
age structures are more filamentary than spherical on the 
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FIG. 4: Configuration dependence of the bispectrum in real and redshift space measured for the dark matter particles in the 
ensemble of LCDM simulations. The four panels show the configuration dependence for fc-space triangles with fe/fci = 2 and 
for scales fei = {0.05, 0.1, 0.5, 1.0} /iMpc -1 . Red and blue colors distinguish between real and redshift space quantities. The 
solid points with error bars show measurements: large solid points are for the real (dots) and redshift space (stars) monopole 
bispectra; the corresponding smaller points are the imaginary bispectra (which should be zero). The solid lines show the 
predictions from HM (thin) and HMs (thick). The triple dot-dash lines show the predictions from PT (thin) and PTs (thick). 
In the top left panel, for clarity we have suppressed the errors on the imaginary bispectra. 



largest scales [3a, |75|. Also, the data appear to be scat- 
tered about the theoretical predictions, with all models 
being equally good fits to the data. 

Considering intermediate scales (k = 0.1 /iMpc -1 ), the 
estimates are much more significant and possess well de- 
fined configuration dependencies - both showing an ex- 
cess of signal for colinear triangles. However, the real 
space bispectrum appears to be in excess of the redshift 
space quantity. This is in contrast to the PT and PTs 



predictions which, whilst qualitatively capture the over- 
all shape, predict the reverse trend, the results being dis- 
crepant by roughly ~ 20%. This problem is also mirrored 
in the HM and HMs predictions, but the natter configu- 
ration dependence of the data is better captured by the 
Halo Model. As was discussed in Sec. IIV CI this flatten- 
ing can be attributed to the impact of the 1- and 2-Halo 
terms acting as effective shot noise contributions. The 
amplitude offsets still require explanation, and we refer 
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FIG. 5: Same as Fig. [4] only this time for the reduced bispectrum. Additionally, where visible the open symbols represent the 
shot noise corrected estimates. 



to our discussion of Sec. I VI Bl for some possible remedies, 
but to that list we may now add the need for loop correc- 
tions to the tree-level bispectra. As was shown in [75| , in 
real space the 1-Loop corrections are significant on these 
large scales. 

On smaller scales (k — 0.5 hMpc^ 1 ), we see that, 
as was noted in Fig. O the real space measurements 
have increased in amplitude and have become much flat- 
ter across the configuration - the HM predictions agree 
rather well with this result. This implies that the statis- 
tic is already in the fully non-linear regime, since the 1- 
and 2-Halo terms are dominating the signal here. There 
is very little relation between standard PT to the mea- 



surements, as expected. Turning to the redshift space 
estimates, it can be seen that the configuration depen- 
dence displays a reasonably strong [/-shape. The PTs 
predictions do not describe this shape very well, but are 
not as discrepant as in real space. However, HMs predic- 
tions capture the form of the configuration dependence 
exceptionally well, but are offset by ~ 10 — 20%. 

Considering the scales associated with the virial radii 
of clusters (fc = 1.0/iMpc ), we see that the estimates 
in real space are surprisingly unchanged and that the HM 
still provides a very good description of the data. For the 
redshift space estimates, we find that there is now a very 
strong [/-shape configuration dependence, in full agree- 
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ment with the results from (35|. Again, the HMs predic- 
tions capture this result remarkably well, although there 
is a small amplitude offset. It is believed that this may be 
mitigated by implementing the improvements discussed 
in Sec. IVTDl and Sec. IVTBl 

We also note that the small discrepancies between the 
real space HM predictions and the simulations, are en- 
tirely consistent with the work of [2(| - in reality haloes 
are triaxial rather than spherical, and this shows up as a 
characteristic increase in signal for colinear triangles and 
a suppression for isosceles configurations. 

F. Code comparison 

Owing to the algorithmic differences between the bis- 
pectrum estimation procedure presented in Section IVI CI 
and that presented in [881 ] . we decided to compare the 
results obtained from these two approaches. Besides 
providing an important cross-check, this also enables us 
to examine how accurately the two codes are recovering 
Q. Overall we found very good agreement between both 
methods, and full results are presented in Appendix [Gl 



VII. CONCLUSIONS 

In this paper, we have provided a new analytic model 
for the redshift space bispectrum of dark matter, haloes 
and galaxies in the plane parallel approximation for the 
redshift space distortion. On large scales, the model pre- 
dictions have a direct correspondence to the non-linear 
perturbation theory and on small scales, the predictions 
are entirely governed by the phase space density of galax- 
ies/dark matter internal to the haloes. This is the first 
time that the information from bulk flows and virial mo- 
tions have been naturally incorporated into an analytic 
model for the higher-order clustering statistics in redshift 
space. 

In our analytic model, the bispectrum is represented 
as a sum over three terms; these correspond to all of the 
possible distinct arrangements of three points in three 
haloes, and we referred to these as the 1-, 2- and 3- 
Halo terms. A practical evaluation of the monopole of 
the bispectrum (a direct observable), with realistic mod- 
els for halo profiles, abundance and clustering, required 
the execution of a set of 4-D numerical integrals. For 
the terms that involved large-scale correlations (2- and 
3-Halo terms) it was shown that these expressions could 
be easily modularized, and so are best computed in paral- 
lel. The 1-Halo term must be integrated with an efficient 
higher-dimensional integrator. 

It was shown that the large-scale predictions in the 
model, which are governed by the 2- and 3-Halo terms, 
can be simplified greatly under the approximation that 
the angle average of the product of anisotropic density 
profiles and large scale clustering of halo centers can be 
performed separately (45l . |46| . This approximation was 



accurate to < 10% on scales k = 0.1/i/Mpc, elsewhere it 
was of the order < 5%. 

The predictions for the bispectrum monopole were 
compared with analytic PT in real and redshift space. 
On very large scales k = 0.05 — 0.1/iMpc , the model 
closely agreed with the redshift space PT predictions, but 
with some small deviations noticeable. It was argued that 
these were due to the 'effective' shot-noise like behavior 
of the 1- and 2-Halo terms in the low-fc limit. On smaller 
scales k = 0.5 — 1.0/iMpc -1 the analytic model showed 
a dramatic departure from the PT predictions and dis- 
played a [/-shaped anisotropy [35l [3a ] . This was the im- 
print in the configuration dependence of the FOG distor- 
tions from non- linear virial motions. No trace of the PT 
remained in the model predictions on these scales. The 
model predictions showed that there was no scale where 
a hierarchical model provided a good description of the 
configuration dependence of the bispectrum. 

The predictions were then confronted with measure- 
ments of the bispectrum and reduced bispectrum from 
an ensemble of numerical simulations. On very large 
scales k = 0.05 /iMpc it was found that the PT and 
Halo Model predictions were equally good, to within the 
errors. On smaller scales, k — O.lftMpc -1 , departures 
between PT and the simulations were noted at the level 
of ~ 10 — 20%. Therefore, studies that use the lowest 
order PT to extract galaxy bias are unlikely to be ro- 
bust on scales k > 0.1/i/Mpc. The Halo Model was a 
better description of the data, but was not in perfect 
agreement. Plausible improvements to the Halo Model 
on these scales were discussed. On even smaller scales, 
k = 0.5 — 1.0 /iMpc -1 , the configuration dependence was 
flat in real space, whereas in redshift space there was 
a very strong [/-shape feature. The numerical results 
were reasonably well reproduced by our halo model pre- 
dictions, a significant improvement over PT that breaks 
down at these scales. 

We cross- validated our results from the numerical sim- 
ulations by comparing our bispectrum estimates with 
those obtained from the independent code of R. Scoc- 
cimarro. The results from the two codes were found to 
be in very good agreement, although our results were a 
factor of 2-3 times more noisy on the largest of scales. 

In future work, we shall extend our analysis to exam- 
ine the halo and galaxy clustering as a function of mass 
and type. It will be also important to resolve whether 
or not the Halo Model's effective shot-noise terms are 
important for modeling real survey data. Owing to the 
fact that different galaxy populations are easily and nat- 
urally included into our description, it is hoped that this 
approach will help to facilitate extraction of information 
from current and future hi-fidelity large scale structure 
surveys of the Universe. 
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APPENDIX A: DUALITY OF CLUSTERING STATISTICS 



This section defines the configuration and Fourier space clustering statistics and their dual relationship with one 
another. It also shows how the cosmological assumptions impose certain important conditions upon these statistics. 

We shall assume that ensemble averages are taken over a volume, V^, of the Universe sufficiently large for the 
fundamental fc-space cell volume to be considered infinitcsimally small; outside of this large volume our stochastic 
fields are exactly zero. This allows us to define Fourier transforms in the volume. We shall also let these fields be 
Ergodic, whence ensemble averages are equivalent to averages over volume. 



1. Real space representation: correlation functions 

The fractional density field of matter is defined: 

*(x,t) = [p(x,t)-p(t)]/p(t), (Al) 

where p(x, t) is the local density at coordinates (x, t) and p(t) is the density of the homogeneous background at time 
t. The n-point connected auto-correlation functions of 6, represent the excess probability from the product of the 
individual independent 1-point distributions of obtaining a particular set of values at all points, e.g. the probability 
of obtaining fluctuations at three points (a, 6, c) can be written: 

P(a, b, c) = P(a)P{b)P{c) [1 + C 2 (a, b) + C 2 {b, c) + C 2 (c, a) + C 3 (a, b, c)} . (A2) 

where C 2 and C3 are the connected 2- and 3-point correlation functions. We may now be clear about what we mean 
by connected correlation function: the connected correlator may not be reduced to sums over products of lower order 
connected correlators. In cosmology these are more commonly written: 

6(x 1 ,x 2 |0 = (5(xi,i)*(x 2 ,t)> c ; (A3) 
&(x 1 ,x 2 ,xa|i) = (8(x u t) ... S(xa,t)) e ; (A4) 
£„(xi,...,x„|i) = (S(xi,t) ... S(x n ,t)) c . (A5) 

These functions obey an integral constraint 

Y J d 3 a;„£„(xi, . . . ,x n ) -> ; n > 1 . (A6) 

This follows from noting that on marginalizing the probability functions over one variable, say the Nth variable, one 
finds that the resulting distribution depends on n — 1-points, and therefore from Eq. (|A2|I must not depend on the 
71-point connected correlation function. 

If the density field obeys the cosmological principle, that is statistical homogeneity and isotropy on scales greater 
than the coherence scale of our fields, then the correlation functions are invariant under translation and rotation of the 
coordinate system. They are also parity invariant real functions and are invariant to exchange of vector arguments. 
Thus: 

£„(xi,...,x„) = £„(xi +x , . . . ,x„ + x ) (Translation); (A7) 
= £„(7£xi, . . . ,72-X^j) (Rotation); (A8) 

= £n(-xi,...,-x„) (Parity) ; (A9) 

= £„(x 2 , xi, . . . , x„) (Exchange) ; 

= £ n (xi,x 2 ,...xi,...,x„) (A10) 
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For anisotropic fields rotation invariance is broken and for inhomogeneous fields translation symmetry is broken. For 
homogeneous fields we may immediately apply the translational invariance and drop one of the vector arguments in 
our function. Setting x = — x„ in Eq. (|A7|) gives, 

£„(xi, . . . ,x„) = £n(xi„, . . .,x (n _ 1)n ) , (All) 

where Xy = Xj— Xj and we shall not write the zero argument in the nth space. In this paper we will mainly be concerned 
with clustering statistics that obey homogeneity, but are anisotropic, as this is exactly the case for the redshift space 
distortion in the plane parallel approximation. Lastly, we have the closure relation: X21 + X32 + • • ■ + xi n = 0. 



2. Fourier space representation: poly-spectra 



Under the conditions stated earlier, the density field <5(x, t) may be equivalently written as an infinite sum over 
plane waves through the Fourier transform, where our Fourier convention is 

d 3 kS(k)e- lkx & <5(k) ' 



(2tt) e 



V,, 



d A x <S(x)e l 



Transforming the density terms in Eqs (|A3IIA5[) . leads to 

P 2 (k 1 ,k 2 )e- ikiri 



£n(ri 



6(ri 2 ) 

^3( r 13) r 23) 
■ • ) r (n-l)n) 



n( d3k 

i=l 



(2tt) s 



[ki + k 2 = 0] 



/ fl{(^} P 3(ki,k 2 ,k3)e- k — * k —3 ; [ kl 



(A12) 

(A13) 
(A14) 
(A15) 

(A16) 

this condition simply arises from the imposed harmonic boundary conditions within our volume: only overlapping 
waves are constructive. The short hand notation [^..J = S (k.% +k 2 + . . . ) has been adopted for the argument of the 
Dirac delta function. The presence of this term ensures that the sum of fc-vectors forms a null vector, ki + • • • + k n = 
and we shall refer to this as the closure condition. For the case of (n — 2), we have the power spectrum P 2 (ki, k 2 ) = 
.P(ki) and for (n — 3) we have the bispectrum P 3 (ki,k 2 ,k 3 ) = P(ki,k 2 ,k 3 ). 

Conversely, the power spectra may also be written as inverse Fourier transforms of the correlation functions: 



i\ l(2^) 3 



P„(ki, . . . ,k„)e" 



-iki -rxn~ 



-ik (n _i)-r (n _ 1)t , 



k 3 = 0] 
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i=l 



where we have very generally defined the n-point spectrum as: 

V;- 1 (5 s (k x ) . . . S s (k n )) = P^( kl , . . . ,k„) [6tn] i^f/v, , 



P(k x ) = J dPrntMe* 1 "" ; 
fl(ki,k 2 ,ks) = [ d 3 r 13 d 3 r 23 t(r 13 ,r 23 )e i ^ r ™ +k * T ™\ 

)n Cn(ri„, • . . , r (n _ 1)n ) e *[ki-n»+-+V C »-x,T (B _ 1)B ] 



Pn (ki , . . . , kn) 



3„ j3„ 



d J ri 



(A17) 
(A18) 
(A19) 



From these relations and the properties of the correlation functions, Eqs (|A7HA10]) , we may now infer the corresponding 
properties for the poly-spectra. Translational invariance in configuration space means that the poly-spectra are 
invariant under a phase shift to the Fourier density fields. Invariance to rotation of the coordinate frame leads to 
rotation invariance of the polyspectra: 

P n (ki, . . . , k n ) = J d\Kv ln ] . . . d 3 [Vr {n _ 1)n ] e„(^r ln , . . . , fcr (n _ 1)n ) e 4^ 1 »+-+kf„- 1) w r(n _ 1)n ] ; 

- f rl\ A 3 *, , t (r- f, } p ,; f r i r " TCTk i + ---+ r ( r "-i)« 7? ' Tk (™-i)l 

= P n {1Z T k l ,...,1l T k ri ) (A20) 
Parity invariance and the reality of the configuration space functions leads to the reality of the poly-spectra: 

P n (k 1) ... ! k„) = P n (-k 1 ,...,-k„) = [P n (k 1 ,...,k„)]* , (A21) 
where the * corresponds to complex conjugation. Many of these properties simplify the analysis in the main text. 
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APPENDIX B: CALCULATIONAL DETAILS 
1. The NFW density profile with the Bullock et al. normalization 

As described in Section IIII Bl to compute the redshift space density profile we require a model for the real space 
density profile p(r) and a model for the 1-point velocity distribution function of particles in a halo. For the density 
profile we adopt the 'NFW model [891 ]: 

p{r) = Pc[y{i + yfY 1 ; y = r/r c . (Bl) 

This model is fully determined by two parameters, p c and r c , a characteristic density and radius. These two parameters 
are not independent, but are related by the mass enclosed: 

_ pA vir c 3 /3 = , R9 ^ 

Pc log(l + c) - c/(l + c) ' " r c ' [tiZ) 

where c is the concentration parameter and is the ratio of the virial radius to the characteristic radius. The virial 
radius is the boundary layer within which all particles have undergone violent non-linear relaxation. It is taken to be 
specified through 

4 

Mvir = -7rr vir A vir/ 9 , . (B3) 

A vir is the density contrast for virialization, which may be estimated from the spherical collapse model. For flat 
universes with a cosmological constant a good fit to the functional form is provided by (90j 

A vir = [(18tt 2 + 82x ~ 39x 2 ] /Q(a) ; x = [Q(a) - 1] . (B4) 

To obtain the concentration parameter as a function of mass we follow the model of Bullock et al. [9l[ . For this 
we have 

a c (M vir ) 

where we take K — 3.0 and where a c (M v [ r ) is the collapse expansion factor for a halo of mass M v ; r . The collapse 
expansion factor for a halo of mass M may be determined through solving the relation 

E±p± (J (FM vir ,a a ) = 5 c , (B6) 
Di(a ) 

where F = 0.001 is taken as a fixed fraction of the initial mass, D(a) is the linear theory growth factor at epoch a. 
The parameter S c — 1.686 is the linearly extrapolated density threshold for collapse from the spherical collapse model, 



where we ignore the slight dependence on cosmology 92]. As was shown by [91( this model provides a very good 



description of the ensemble average properties of dark matter haloes. More sophisticated models may be constructed 
that take into account that haloes are more complicated, e.g. a halo of mass M drawn at random from the ensemble 
will have a concentration parameter that is drawn from a probability distribution of possible concentrations. In 
addition one may include sub-structure [§3| or halo triaxiality [2(| UM- However, we shall leave these additional 
embellishments for future study 

2. Conversion between Sheth &; Tormen and Bullock et al mass deflntions 

The definitions of halo mass that are used in the Sheth & Tormen mass function and the Bullock et al. model for 
the density profile are inconsistent. We recall that the Sheth & Tormen halo mass is defined: 

4 
— ; 

3 



M ST = -7rr| T 200p . (B7) 



We resolve this inconsistency using the methodology of [72j, which briefly is as follows: For a given halo of the NFW 
type, the physical values of the characteristic density and radius are independent of our specific choice of halo mass. 
Using the relation for the physical density as a constant we arrive at the mapping 



CST 



200 



log(l + c vir ) - c vir /(l + c vir ) 



log(l + cst) - cst/ (1 + cst) 



(B8) 
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Thus if we take a Bullock et al. mass and derive the appropriate c v i r we may then solve the above expression to find 
cst- Following this we may then obtain the corresponding Sheth & Tormen mass through application of the relation, 

200 /c ST x "' 



M ST = -7— — M vil . (B9) 



3. 1-Point velocity distribution profile 



For the ID velocity distribution function we adopt the standard Maxwellian distribution [45|, [46J, l47|, \9_ 



V[u z |ai D (M vir )]dM 2 = 



27TCT 



■ exp 



ID 



2of D 



du z 



(BIO) 



where oid is the ID velocity dispersion. For haloes that possess an isothermal density distribution, this quantity is 
related to the halo circular velocity (V c ) through the following relation [95j 



'ID 



(Mvfr) = eV c 2 /2 ; 



GM V 



(Bll) 



Note that we have included a parameter e into Eq. (|B11[) , this may be used to account for the fact that the relation is 
only approximately true for the NFW density profile model. It also serves the further purpose of allowing us to turn 
off the fingers-of-god through setting e — ► r] -C 1. As discussed in Section IVI Bl and shown in Fig[5]e = 0.76 provides 
a reasonable fit to the velocity dispersion mass relation from simulations. On combining the above relations we have 



'ID 



(M vir ) 



100 e 



n m (o) H(a) 

' V 



(B12) 



Notice that the ratio <7iD(-M v i r )/7Vi r is independent of halo mass. One immediate consequence of this is that in 
redshift space the ratio of the line-of-sight projection of particles in a halo compared to the transverse length will be 
a constant, am / (Hr v i T ) ~ 3, regardless of mass. In other words, FOGs lead to density profiles which are self-similar. 
In our analysis we take e = 0.76 (See Fig. [5]). Finally, the Fourier transform of the velocity distribution is 



V(/iifci 



exp 



1 



[ki[ii<7i D (M vil )] : 



(B13) 



APPENDIX C: EULERIAN PT AND HALO-PT KERNELS 



The first two symmetrized Eulerian PT kernels for the density and divergence of the velocity field are 

G\.i 



5 M12 

7 2 

3 M12 

7 2 



Si + Si 

92 qi 
Sl + Sl 

32 qi 



2(Mi2) 2 
7 

4(Mi2) 2 



7 



(CI) 
(C2) 



where Fi,...j = -F)(qi, ■ ■ ■ , cji) and where = qi • q 2 /<7i<?2- 

The Halo-PT kernels, symmetrized in all of their arguments, may be written in terms of the Eulerian PT kernels 
up to 2nd order as [38(: 



F, 



he 



b (M) ; 



Fr = b 1 (M)W(\k\R)F 1 ; 
Fi c 2 = h(M)W(\k\R)F 1:2 



b 2 (M) 



WQqxlfyWQqaWFiFa 



(C3) 
(C4) 

(C5) 



where F^ c • = c (qi , . . . , q^ \M, R) and where k = qi 
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10 14 10 15 
Mass [M Q /h] 

FIG. 6: Mass verses ID velocity dispersion measured for FoF haloes in simulations. Points show mean and 1-sigma errors for 
measurements from the numerical simulations. The solid and dash line shows the predictions from Eq. lB12l with e = {0.76, 1.0}. 
The bottom panel shows the ratio of the data with respect to the e = 0.76 model. 



APPENDIX D: ROTATION MATRIX 



Owing to there being several equivalent ways to define the Euler angles and hence the rotation matrix TZ(ji, 72, 73), 
we make explicit our adopted choice. The angles are defined as follows: 71 describes a rotation of the coordinate 
system around the z-axis; 72 a rotation around the new y'-axis; and 73 a rotation around the new z"-axis ([96]). 
Thus, the components of any vector k specified in some initial Cartesian system can be transformed into the scalar 
components of the new rotated basis vectors through k' = 7^.(71,72, 73 )k. The z — y' — z" rotation matrix is (9d| : 



([C72 C71C73 - 571 S73] [C72 571 C73 + C71 573] -57 2 C7 3 \ 
[-C72 C71 S73 - 571 C73] [-C72 571 S73 + C71 C73] Sj 2 Sj 3 , (Dl) 
S"72 C71 6*72 Sji C72 J 



and we employed the economic notation Cx — cos x and Sx = sin a;. 
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APPENDIX E: REDSHIFT SPACE POWER SPECTRUM MONOPOLE IN THE HALO MODEL 



The redshift space power spectrum of tracer particles a, can be written in the linear halo model as 

^(k) = ^, 1H (k) + ^ )2H (k) : 
^,m(k|i?) 1 



^iH(k) 
^, 2H (k) 



[W{k\R)Y 

^, 2H (W _ i 
[W{k\R)f Pi 



4 [ dMn(M)[W a ] 2 \U s (k\M)\ 2 , 

Pa J 



i=i 



P{ c (k\M 1 ,M 2 ,R) 
[W(k\R)f 



(El) 
(E2) 

(E3) 



At linear order the redshift space power spectrum of halo seeds is: 
P£ c (k\M u M 2 ,R) = Z 1 (k\M 1 ,R)Z 1 (k\M 1 ,R)P 11 (k) 



[^(fci?)] 2 6 1 (M 1 )6 1 (M 2 )P 11 (fc) {1 + [Pi + M + /?i/3 2 /i 4 } ; Pi 



m) 



(E4) 



where 6 — can be seen from the fact that the halo and density fluctuation fields are by definition mean zero 
fields and recalling that at linear order ^<5{ l (Y|M)) = bo(M) + b\(M) (5\{r\M)). The redshift space power spectrum 
monopolc is thus 

1 



KM k ) = ^ I dMn(M)[W a r\U a (k\M)\ 2 n ( ^[ka(M)} 

Pa J 

i r 2 

P«M k ) = — ]]_{dMin(M i )bi(Mi)[W a }iU(k\M i )}n^{ka 2 (M)] 

Pa J 



where we have defined the redshift space multipole factors 



21 + 1 
2 

21 + 1 



d[iVi{(i) exp [-a 2 v 2 ] ; 
J diiVi (m) [1 + A^i 2 + B^} exp [-6 V] 



(E5) 
(E6) 

(E7) 
(E8) 



Here a 2 = nk 2 a 2 (M)/2, ; b 2 = k 2 [a 2 (Mi) H h<r 2 (A/„)] /2 we have set A = + fa a nd B = fofo, and we 

have assumed our Gaussian model for the 1-pt velocity distribution function from Eq. (|B13|) . Thus, the monopole 
[I = 0; Vo(fi) — 1] moments are 

erf [a] 



2 a 



,Jb) - 7 ^l^ + 2b 2 A + 3B]-^^l 2i r [ A --B l --,:U^ 



K ( 2 °Ub) 



(E9) 
(E10) 



Our expression differs from that of [45l . l4q | , but is consistent with the formulation of [47j . For further discussion and 
comments on this subject, and for an evaluation of the power spectrum to higher order in the Halo-PT series, see [4^ |. 
Note that when a <C 1 and b <C 1, then 



K { °i(a) = exp(-a 2 ) V , V - a 2 ^ = V ,) ±r , a 



^ (27 H- 1)1 



K 2 %{b) 



exp(— 6 2 ) 
464 

exp(— 6 2 ) 
464 



1 



[46 4 + 2b 2 A + 35] ^ 2 ' + !)!! ^ ~ 2b ^ A + B ) ~ 3S 



(Ell) 
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[46 4 + 26 2 A + 3B] 



2b 2 46 4 



15 105 
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1 H 
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46 4 
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5" 




IT 




35 


+ 15 


7^ 
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♦-) 



3 5 

reducing to the Kaiser formula on large scales [4C 



(E12) 



25 



APPENDIX F: IMPACT OF SHOT NOISE ON THE REDUCED BISPECTRUM 



It is of interest to consider how standard shot noise and also the halo model effective shot-noise terms impact the 
reduced bispectrum. On large scales Q HM can be written: 



B 



PT 



iHM 



^ [Pi + P-2 + P 3 ] + -J- 

2H,B "lH.E 



P 



«1H,P 



2cyc 



(Fl) 



where we have added a sub-script P or B to distinguish between shot noise terms from the halo model power spectrum 
and bispectrum, respectively. 

In the low sampling limit n^P <gC 1, we have the result that 



Q 



HM 



{[W a f)([W a ]) 

(\w a ?) 2 



(F2) 



for the case of standard shot noise, the term in square brackets is unity and we have Q A = 1/3, where super-script d 
means discrete. 

In the high sampling limit, n^P ^S> 1, the denominator in Eq. [FT] becomes, 



1 



; fac 



2 P1+P2+P3 



1 



Q 



fac 



P Qfac 



(F3) 



where Qf ac = P1P2 + 2cyc and where we have treated the last two terms in the square brackets as small quantities, 
relative to Qi ac . On replacing this in Eq. (|F1[). we find 



Q 



HM 



Q PT = 



Pi 



Ps 



Q 



fac 



1 



-(i + Q PT ) 



1 



/fac 



1 



'1H,P 



-(i + Q 



PT\ 



(F4) 



We can also use the above expression to derive the effect of standard shot noise on the reduced bispectrum, by simply 
considering all of the number density terms to be identical, and this gives 



Q d -Q 



PT 



(Pi + P 2 + P 3 ) 



nQ 



fac 



(1 + 2Q 



PT\ 



1 



n 2 Q 



fac 



-(2 + 3Q 



PT\ 



(F5) 



Considering standard shot noise first, in Eq. (|F5|) . we see that the effect of discretization of matter is always to reduce 
the value of Q d relative to Q PT (the continuum limit case) . Considering now the Halo Model, in Eq. (|F4[) , we see that 
the difference between this and Q PT depends on the sign of the quantities in square brackets. For dark matter, the 
first term can be seen to be negative, since ?Z2H,b = %h,p- However, the sign of the second term is not as obvious to 
deduce. If it is negative, then the effect is as for standard shot-noise; on the other hand, if the reverse is true, then 
Q d > Q PT . 

Lastly, the configuration dependence of Q d — Q PT in the standard shot noise case can be understood from the 
following: if we assume that all fc-vectors are larger than the turn-over scale in the power spectrum, then the quantity 
Qfac(^i2 = 1) > Qfac(#i2 = 0). This implies that the difference is largest when k\ and &2 are parallel. 



APPENDIX G: BISPECTRUM CODE COMPARISON 



In this appendix we compare results from our recipe for estimating the reduced bispectrum, presented in Sec. IVI CI 
with those obtained from an independent prescription used by one of us over the years, e.g. [88| . which uses full 
sampling of all fc-space triangles on the Fourier grid. The two methods are very similar, but some subtle differences 
exist, these can be summarized: for the "full sampling code": 

(1) B is estimated in linear bins of thickness a 1 X kf for k — 0.05, 0.1[/iMpc -1 ] and 4 x fc/ for the case k = 
0.5, 1.0[/iMpc _1 ] [kf = 2-7T /L); whereas for our code estimates are made in Alog 10 k bins of thickness 0.05; 
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FIG. 7: Comparison of the reduced bispectrum estimates obtained from the method presented in Sec. lVI Cl with those obtained 
from the "full sampling code" of |88| . As in Fig. [S] our results are represented by solid symbols, and those from the alternate 
method are denoted by open symbols. Again, errors are on the mean and real and redshift space estimates from the same code 
have been slightly offset in the :r-axis to enhance clarity. 



(2) the configuration dependence of B is estimated as a linear function of ks, whereas for our code it is estimated 
as a linear function of #12; 



(3) all of the available independent fc-modes are used, whereas we sub- or over-sample modes from the available set 
depending on the number of available modes; 



(4) Q is constructed from estimates of B and P in a post-processing fashion, whereas we estimate Q on-the-fly for 
each triangle that is used in the estimate. 



Fig. [7] shows the results of this comparison. The solid symbols denote our results, and the corresponding open 
symbols denote the results from the "full sampling code". Overall we find very good agreement between both 
methods. On the largest scales that we have considered, k\ — 0.05 ft,Mpc _1 , it appears that our method is a factor 
of 2-3 times more noisy than that of R. Scoccimarro, however we have a factor of 2 more bins in #12, which accounts 
for some of this discrepancy. On smaller scales ki > 0.1 /iMpc -1 the estimates are of comparable quality, with ours 
being slightly more noisy. The discrepancy on large scales owes to the fact that we have sub-sampled triangles from 
the possible set, this can be mitigated by oversampling from the number of available modes. On smaller scales the 
benefits of our approach are that we may obtain a high accuracy estimate without requiring all of the triangles and 
this also has the practical advantage of keeping the computational time tolerably low. 
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